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Abstract 

We describe effects of anisotropy caused by the crystal lattice in d-wave 
superconductors with s-wave mixing using the effective free energy approach. 
Only the d-wave order parameter field d is introduced, while the effect of 
the s wave mixing as well as other effects breaking the rotational symmetry 
down to the fourfold symmetry of the crystal are represented by a single four 
(covariant) derivative term: 7]d*(Ily — li^)'^d. The single vortex solution in a 
phenomenologically interesting range of parameters is almost identical to the 
two order parameters approach. We analytically consider the most general 
oblique lattice and orientation, but find that only rectangular body centered 
lattices are realized. A critical value rjc at which a phase transition from the 
rectangular lattice to the square lattice takes place. The influence on the phase 
transition line is discussed. The formalism is extended to the time dependent 
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anisotropic Ginzburg-Landau equations in order to calculate the effect of the 
anisotropy on the flux flow. The moving vortex structure is established and 
the magnetization as function of the current is calculated. Although the linear 
conductivity tensor is rotationally invariant due to the fourfold symmetry, 
the nonlinear one shows anisotropy. We calculate dependence of both direct 
and Hall I-V curves on the angle between the current and the crystal lattice 
orientation. 



I. INTRODUCTION AND SUMMARY 

It is widely believed that the superconductivity in layered high Tc cuprates is largely 
due to the d(^x2_y2-^ pairing The evidence for the d - wave pairing comes from variety 
of different measurements. A partial list includes the /xSR measurements of the penetra- 
tion depth 0, quantum phase interference 0, angular resolved photoemission thermal 
conductivity 0, vortex lattice structure observed using neutron scattering ^ and tunnel- 
ing spectroscopy ||^ and nuclear spin relaxation rate |p. While most of these experiments 
directly probe the energy gap and the low lying excitations, the vortex lattice observations 
are different. They try to relate the structure and interactions of Abrikosov vortices to the 
nature of the order parameter. This would be extremely difficult if the order parameter 
would be a pure d wave. The structure of a single vortex and even the vortex lattice on the 
level of the "macroscopic" Ginzburg - Landau equations would be the same as for that of 
the usual s wave superconductors. 



There are numerous indications, both theoretical [T^M and experimental [plp[, that even 



though the major pairing mechanism is the d wave pairing, there is a small admixture of 



the s wave pairs in the condensate. Ren et al ||rT| using a phenomenological microscopic 
model (in the weak-coupling limit) and Soininen et al [|I2| considering attractive nearest 
neighbors interaction derived a two field effective Ginzburg-Landau (GL) type theory. The 
two complex order parameter fields, s and d describe the gap function in corresponding 
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channels. The most general free energy of this kind reads |T3|,|T4[: 

+7s|ns|2 + 7,|nrf|2 + 7.[^*(n^ - K)d + c.c] (i) 

where 11 = — iV — e*A is the covariant derivative and e* is the charge of the Cooper pair 
(throughout this paper we use the convention c = h = 1). Within a particular microscopic 
model there might be some relations between these coefficients, but since the ultimate mi- 
croscopic theory is not known as yet, all of them should be considered as phenomenologically 
fixed parameters. 

Using equations following from this free energy or more fundamental equations (see recent 
quasiclassical Eilenberger equations treatment in one obtains a characteristic four- lobe 



structure with four zeros for the s - wave inside a single vortex ||13|JT6[1 . Therefore the vortex 
core looses the full rotational symmetry and only the fourfold D4,h symmetry remains. The 
distribution of the magnetic field was also obtained recently [|I6],|13|. A somewhat different 



structure, however, was found in another recent numerical work |[TB| (our analytic calculation 
confirms that of [|TB],|T^ and contradicts that of [0). 



Outside the core the s-wave vanishes, while the d-wave becomes rotationally invariant, 
indistinguishable from the usual Abrikosov solution. Therefore, to look for differences in 
the behavior of vortices one would like to be closer to ifc2; so that the core will be more 
important (which is not easy for high superconductors due to their large Hf.2)- However, 
since the fourfold vortex core structure comes into a conflict with the high symmetry of the 
triangular lattice, the asymmetry of vortices can distort the usual triangular vortex lattice 
at already accessible flelds much lower than Hc2- Another phenomenon in which the vortex 
core plays a major role is the dissipation in the course of the flux flow. These are the two 
major phenomena in clean superconductors in which one can look for anisotropy effects (we 
neglect pinning and other disorder and fluctuations, and concentrate only on YBCO single 
crystals). 

In this paper we study in detail the above two phenomena, vortex and the vortex lattice 



structure and the flux flow, using a greatly simplified model: time independent and time 
dependent one component effective Ginzburg - Landau equations. This formulation, essen- 
tially without loss of generality, allows us to avoid numerical methods and to extend and 
clarify various delicate questions about the single vortex and the vortex lattice structure for 
which there is some controversy or lack of concrete proof. In particular the location of the 
phase transition to the square lattice is calculated and depends just on one combination of 
the parameters. Moreover the quantitative discussion can be extended to study moving flux 
lattices, which, as is well known [l^JS^ , are much more demanding, as far as the calculational 



complexity is concerned. Their shape and orientation of moving lattices (or large bundles) 
with respect to the crystal lattice and electric fleld (or current) is determined. Then we 
calculate the current due to the flux flow as function of the external electric fleld. Since the 
deviations from the full rotational symmetry in the flux flow can be clearly seen only in the 
nonlinear regime (the fourfold symmetry forces the full rotational invariance of the Ohmic 
conductivity tensor) one has to go beyond linear response.. The simplifled formulation is 
indispensable in this case |T^J2n], but the result turns out to be remarkably simple. In the 



rest of this section we outline the basic assumptions, methods and results pointing out where 
in the paper more details can be found. 

Within the two fleld formulation Soininen at al [l^ observed that in predominantly 
d - wave superconductor the s - wave component is generally very small: it is "induced" 
by variations of the larger d component. Mathematically the dominance of the d - wave 
follows from the fact that coefficient of the d*d term, —ad = a'{T — T^), is negative , while 
that of the s*s, a^, is positive in the free energy, Eq.(||). Therefore, it is the d fleld which 
acquires a nonzero value. Then the rotationally noninvariant derivative term s*(n^ — n^)(i 
" communicates" the deviations from the condensate value of d to s. Note that this is the 
only term up to (scaling) dimension three in the free energy which breaks the full rotational 
symmetry. If its coefficient 7^, is not very large, the s fleld never becomes comparable to d. 
Soininen et al |]TB|| observed that even near the core, where d is the smallest, it is nevertheless 



larger then s by a factor of 20 at least. This in particular means that many small terms 



like |s|^are irrelevant. The field d is the critical field near T^, while s is not and can be 
"integrated out" perturbatively. It will be explained in some detail in Section II, that this 
generates an effective (scaling) dimension five term for d, breaking the rotational symmetry, 
so that our starting point to study the rotation symmetry breaking effects is an effective 
free energy 

Here we have replaced 7^ by a more conventional notation l/2md and parameter rj = '^l/as 
quantifies the deviation from the rotational symmetry. Its relation to other parameters in 
the two field approach is derived in Section II. We calculate all the above mentioned rota- 
tional symmetry breaking effects to the first order in 77. This formulation follows the general 
philosophy of effective free energy written in terms of critical fields only. The noncritical 
fields just renormalize the coefficients. The rotational symmetry breaking term has dimen- 
sion five, but breaks the symmetry and is therefore a "dangerous irrelevant term" using the 
terminology of critical phenomena [^. The contributions to it might come not only from 
the d - s mixing, which always give positive 77, but also from other sources. Even in conven- 
tional superconductors such effects exist [^. In YBCO, twinning might be an important 



contribution to it. This formulation avoids the problem of the second phase transition at Tg 
that one encounters assuming = a'{Ts — T) in the two field formalism. 

The single vortex solution is obtained in Section III. It is almost identical to the solutions 
obtained earlier within the framework of the two order parameters theory. One still can define 
the "effective s-wave field by s = —^(11^ — 11^) (i and observe the four lobe structure, see 
Fig. 2 and 3 for d and s components respectively. Relation to earlier work (discrepancies or 
common points) are summarized in the Appendix. 

The vortex lattice near Hc2 is studied comprehensively in Section IV. The simplicity of 
the formulation allows for an analytic study of all the possibilities, not considered before or 
considered using uncontrollable approximations. The degrees of freedom we include in the 
analysis contain: (1) an arbitrary rotation angle if between the crystal lattice and the vortex 



lattice (Fig. 4), and (2) all the possible lattice, not only the rectangular ones considered 
before ( [|1^], |T^, |18|). Instead of using the variational method to solve the linearized 
set of equations, we solve it exactly even in the moving lattice case. This is the first time 
that the lattice is demonstrated to be rectangular body centered using the most general 
lattice in the analysis. We tabulate the lattice characteristics for different rj in Fig. 7. At 
certain value of rj there is a phase transition from rectangular to a more symmetric square 
lattice first noticed [|l3l . The exience of a phase transion becomes obvious in our formulation 



in which the effective strength of the four fold symmetry is proportional to the magnetic 
field, characterized by a dimensionless parameter rj' = rjmde*!!. In low fields, the four fold 
symmetry is subdominant, so the lattice is closer to the triangle lattice. In high fields, the 
four fold symmetry dominates, so the lattice becomes square. We find that the transition 
occurs at t]'^ = .0235. The upper critical field line Hc2{T) is given in Eq. p7|) . However at 
the end of this Subsection IV. C we caution against too direct interpretation of this result by 
considering other possible effects which are isotropic and also contribute to the curvature. 

The moving lattice solutions are derived in Section V from a time dependent generaliza- 
tion of Eq.(|^). They are not only needed for the nonlinear conductivity calculation which 
follows, but are also interesting in their own right, since they are, in principle, observable. 
In the one field formalism there is only one additional parameter to be added to take the 
dissipation effects into account: the coefficient of the time derivative term of d. Unlike in 
the case of the pure s-wave superconductor, the moving (with arbitrary, not infinitesimal, 
velocity) lattice solution in this case can not be obtained from the static one by a simple 



Galilean boost pO[. It is a nontrivial problem and we were able to solve it perturbatively in 
1]. Unlike the s - wave moving lattices (which are triangular [0), orientation is determined 
by the direction of the crystal lattice as well as by the current direction. The "dynamical 
phase transition line as function of current and its orientation with respect to the atomic lat- 
tice are quantitatively discussed in Subsection V.B and the result is given in Eqs.(|6^, |6^ , |6^ ). 
Magnetization close to Hc2 (or the Abrikosov /?^) is a very simple function of the current J 
and its direction. The result for the Abrikosov Pa is: 
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/3a (E) = + ajE^ I cos 29 1 , (3) 

where is an angle between the electric field E and an axis of the underlying atomic lattice 
axis. The number f3^ (typically a bit larger then 1) is the usual Abrikosov (3 parameter for a 
given lattice defined in Eq.(^T]) and the constant c is given in Subsection V.C (Eqs. (|75| , |76D ). 

Corrections to the linear conductivity tensor (which cannot break the rotational sym- 
metry) are briefiy discussed and a detailed calculation of the effect of the anisotropy on the 
nonlinear fiux fiow are given in Section VI. The result is remarkably simple. In addition to 
isotropic linear part there is an anisotropic cubic in E term is: 

AJ = ^^^(l + cos4e) (4) 

and the Hall current is 

^Jnaii = -V ^Og.ff4 sin 49 (5) 

The two nonlinear contributions to currents are simply related. In these expressions 7 is 
the coefficient of the time derivative term in time dependent Ginzburg - Landau equation 
Eq. (p!OD . Both direct and Hall I-V curves depend on the angle between the current and 
the crystal lattice orientation via the fourth harmonic only. The result, contains only cubic 
dependance of the currents on the electric field, higher orders being cancelled. 

Finally in Section VII we conclude by briefiy discussing possible experiments to observe 
various above mentioned effects, as well as some generalizations. 

II. TIME INDEPENDENT AND TIME DEPENDENT EFFECTIVE GINZBURG - 

LANDAU EQUATIONS. 

On very general grounds superconductivity is sometimes considered as a phenomenon of 
"spontaneous U{1) gauge symmetry breaking " irrespective of the mechanism of pairing 
or channels in which it occurs. The U{1) electric charge symmetry should be represented 
by a single order parameter: the superconducting f/(l) phase. Mechanisms of superconduc- 
tivity, as far as connection to the gap functions appearing in the microscopic description 



is concerned, may differ, but this general order parameter representation of the supercon- 
ducting phase remains the same. While in pure s - wave or d - wave superconductors the 
f/(l) phase is simply identified with the phase of the gap function, in more complicated 
microscopic theories with few channels opened, the superconducting phase is just the com- 
mon phase of various gap functions. Therefore, quantities other than the phase which enter 
various phenomenological Ginzburg-Landau (GL) type equations, although useful, are not 
directly related to the spontaneous gauge symmetry breaking. 

The s-d mixing two component GL free energy Eq.(||) leads to the following set of equa- 
tions: 

(7^0^ -ad)d + 7,(nJ - Ul)s + 2(32\d\^d + (33\s\^d + 2/?4s'rf* = (6) 
(7,n2 + as)s + 7,(n2 - Ul)d + 2(3i\s\h + f33\d\^s + 2Pidh* = (7) 

As we discussed in the Introduction, the noncritical s component is induced by the 
d — s mixing gradient term. Therefore the length scale of the variations of the field s is 
the d = wave's coherence length = ./^. Consequently the derivative term in Eq.flTD 
7sn^s ~ ils/^d) ^ is small compared to agS. This requirement <^ 1 (typically ^ ~ 1) 
holds for the vortex solution of ||T3| , p!7|Jl6[| and is, in fact, an excellent approximation in both 
near and far from the core regions. Then, according to Eq.(|^), the field s is, to the first 
order in l/a^: 

s^-^{Ul-Ul)d (8) 
as 

Substituting this equation back to Eq.@, we obtain, to first order in l/cts, 

1 _0 \ . O -r2\2 



W -ad)d-7]{Wy-Uiyd + 2(3\d\'d = (9) 



.2md 

where 7^ was replaced by a more standard notation and rj = The second term should 
be treated as a perturbation. In principle there are terms of higher order in l/a^, but one 
cannot take them consistently into account without simultaneously including additional 
terms in the original two field GL equations. From Eq.(|]) the effective free energy Eq.@ 
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follows. Of course the effective free energy f^ff Eq.(|^) is valid only if possible higher orders 
of the field d and derivatives are neglected. 

Note that even the linearized set of Eq.(|^,0) is highly nontrivial. Authors of Ref. |13 



resort to the variational estimate to find a solution. On the other hand, the linearized Eq. 
can be easily solved perturbatively in l/a^. Another advantage of this equation, especially 
as far as relation to experiments is concerned, is that the number of coefficient is much 
smaller. Instead of 6 additional parameters in the two field free energy, there is just one 
additional adjustable parameter compared to the usual s - wave GL equations. 

The effective free energy approach can be also motivated by considering the fluctuations 
of the two fields theory. Assuming that only the d field is critical, one can integrate out the 
s field fluctuations perturbatively. The lower order terms coming out from this analysis are 
precisely the same as /g// as expected. In principle, the coefficients of the effective one field 
free energy should be fixed by a microscopic theory in the same spirit as the way that the 
coefficients of the two component equations should be fixed. The general form of the effective 
free energy can be obtained just by the dimensional analysis and symmetry. Generally, we 
should allow all the terms invariant under the group -D4/1 with dimensions five or less. It is a 
well known property of the 1^4/4 symmetry that at the level of the dimension three (relevant) 
terms, full rotational symmetry is restored. 

Therefore, one has to consider "irrelevant" (scaling) dimension five gradient terms in 
order to break the rotational symmetry down to D^^. Apriori there are five such terms: 
d*Il'^d, d*Ilyd, d*lllllyd, d*Il^Ilyd, and d*lllllyd. The last two break the reflection sym- 
metry, while out of the remaining three, one can only construct two linear combinations 
invariant under the rotation of it/2. They are d*{Ily — Iliyd and d*{lPyd. The second term 
is fully rotational invariant, and therefore is not important for studying anisotropy effects. 
It, however, contributes to effects not directly related to anisotropy such as the shape of the 
phase transition line. We will come back to this point in Section IV. Similarly, all the other 
dimension five terms that one should in principle include are rotationally invariant. They 
just add a small contribution to the rotationally invariant parts of physical quantities and 



need not be included when studying the rotational symmetry breaking effects. 

The time dependent GL equation in the one field formulation is also quite simple: 

7 + ^e*$j d=- (^n^ -a,^d + r]{Ul - UD'd - 2(3\d\% (10) 

where $ is the electric potential. The above equation describes the time evolution of the 
order parameter and will be used to describe moving vortex systems. It involves only one 
additional parameter 7 compared to the 2x2 matrix for the two field formalism. Although, 
in principle, this parameter describing various dissipation effects can have a complex part 
|[T9| , we will consider only real values. 

III. THE SINGLE VORTEX SOLUTION 

In this section we find an isolated vortex solution of the one component equation Eq. 
near Hd-The opposite case of magnetic field close to ifc2 will be considered in next section. 
We measure the field in units of the vacuum expectation value \E'o = \fc(d/2f32, and length 
in units of the coherence length = l/y/2'mdad ■ In strongly type II materials (when the 
Ginzburg-Landau parameter k is very large), as is the case in high superconductors, we 
can safely ignore the magnetic field and the dimensionless GL equation becomes: 

(-V' - l)d - A(V,J - Vlfd + \d\^d = 

where A = Arjm'^ad is our dimensionless small perturbative parameter.. We solve it per- 
turbatively in A as follows. Let d = do + Xdi , where do = /o(r)e*'^ is the solution of the 
standard unperturbed GL equation. Then the first order equation in A is 

(-V' - l)di + {2\do\^di + dldl) = (Vj - Vlfdo (11) 

The angular dependence of di is easily observed to contain only three harmonics: e"^**^, e'^^'^ 
and e^*"^. This is consistent with the fourfold symmetry which is built into the theory. We 
therefore decompose di into combination of these three harmonics: 



10 



(12) 



The equation becomes 



'A 

d 



1 d 
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+ ^ + 1 /-3 - /o (2/_3 + /s 

Id 25 



1 d 

dr"^ ' r dr 



+ -^-12+1/1-3/0/1 



-^-3(r) 

-^5(r) 

-Ji(r) 



(13) 
(14) 
(15) 



with Jj defined by: 



(V; 



72\2 



/o(r)e^ 



e*Vi(r) + e-3'<^J_3(r) + e'^*V5(r 



As is well known, the analytic expression for /o does not exist, however, there are a number 

, the set of linear 



of known good approximations. Using one of them /o(r 



equations are then solved numerically (the third equation decouples from the first two). The 
results are shown in Fig.l. The d - wave configuration is basically indistinguishable from 
that of the two fields formalism for A = 0.15, see Fig. 2. 

Note also that within the same approximation and normalization, the s component is: 



. = A'(V^-V^)cio 



(16) 



where A' = 2'^^md{ad/ as) 



A 



2-yy rrid 



is another dimensionless small parameter. It's easy to 



see that s has the asymptotic behavior 



s ~ re' 



-i<f> 



(17) 



The s field is plotted in Fig.3. The different winding numbers in near and far asymptotic 
regions give rise to four additional poles in the s component in the intermediate region. This 



confirms calculations of ||T3[ even though some asymptotic analytic expression used there 
to obtain the numerical results disagree with ours. The comparison with [|1^] and [jl^ is 
presented in Appendix. 



11 



IV. THE VORTEX LATTICE NEAR Hc2 



In this Section we follow a generalization of the Abrikosov's procedure [0,^ to inves- 
tigate the structure of the vortex lattice near Hc2- One first ignores the non-linear terms 
in the GL equation and finds a set of the lowest energy solutions \l/fc„(a;, ?/)of the linearized 
equation. The lattice solution is constructed as a linear superposition 

d{x,y) = Y.Cn^kA^,y) (18) 

n 

in such a way that it is invariant under the corresponding symmetry group of the lattice. It 
is well known that the free energy near Hc2 is monotonic in the Abrikosov's parameter (3a, 
which is defined by 13a = In the last step needed for some applications, the 

overall normalization of the order parameter is variationally fixed by minimizing the free 
energy including non-linear terms. 

A general lattice in 2D can be specified by three parameters a, h and a, where a and 
h are the two lattice constants, while a is the angle between the two primitive lattice vec- 
tors (Fig 4). Flux quantization constrains them, so that Habsma = $0- In the d-wave 
superconductors the rotational symmetry is broken, therefore the relative orientation of the 
vortex lattice to the underlying lattice becomes important. Later we will denote (p to be the 
angle between a and one of the axes of the underlying lattice. In Abrikosov's original paper 
01 he assumed C„ = Cn+i and obtained the square lattice, later Kleiner et al generalized 
the procedure to the case where C„ = Cn+2- In this way all the rectangular body centered 
lattices can be included in the analysis. In previous works on d-wave superconductivity ref. 



13,17], the same formalism was used, however, this does not include the most general 



lattice. In this section we follow a more generalized formulation of Ref. p5[ which covers 
all possible lattice types. 

In our case we first solve the linearized GL equation pertubatively in the anisotropy 
parameter rj. And then we obtain an analytic expression of (3a as a function of a, h (or a) 
and if . Finally, we minimize the free energy analytically with respect to if and numerically 
with respect to a, h (or a) to find the lattice structure. 
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A. The perturbative solution of the Unearized GL equations 



We start from the effective linearized GL equation Eq. 

1 



U'd - (n; - Uiyd = aad, (19) 



2md 

where for later convenience we have moved add to the right hand side. It is important to note 
that in Eq. (^) we have assumed that the coordinate system and the underlying microscopic 
lattice coincide. Later it will be convenient to orient the coordinate system {x, y) with the 
Abrikosov vortex lattice rather the atomic crystal. In general, if the crystal is rotated by an 
angle ip counterclockwise with respect to the coordinate system, Eq.(|19|) becomes 



^U^d - 7] [cos 2^ {III - nj) + sin 2^ {lijiy + n^^n^)] ^ d = add. (20) 
It is convenient to introduce dimensionless creation and annihilation operators defined by 



H, 



= — 
V2 



V2 

'h, 



where Il± = 11^; ± illj^ and the scaling parameter Ih = l/ve*iJ is the magnetic length. In 
terms of a and , Eq. ([T9D becomes 



ata + - - V (e+2^^at2 + e-^'^a^ 
2 



dix,y) = —d{x,y). (21) 



Here dimensionless parameter rj' is given by t]' = rjmde*!!. For later convenience, we have 
defined an unperturbed (conventional) upper critical fields = $o/(27r^J) = 2mdad/e*. 

It is convenient to choose the Landau gauge A =Hxy, since in this gauge the momentum 
operator p^^ commutes with d. Therefore we can choose d{x, y) to be an eigenvector of Py 
with eigenvalue k: d{x,y) = exp{iky)tljk{x). The operators a and in this gauge become 



/ d 
^[2 \dx 



« = 4(:7- + 5^h (22) 
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where x = (x — xq)/Ih is dimensionless with xq = klj^. Standard perturbation theory for 
the Schrodinger type equations gives for the lowest eigenvalue: 

§^ = l-2v' + 0{r^'), (24) 

This determines the upper critical field. Note that the relative angle (p does not affect H 
in the lowest order. We will come back to examine this result more closely later. The 
corresponding eigenfunction iIj{x) is: 

V'(5) =^Po + V'i^i = |0) + v'-r"^- |4) + Oiv'') (25) 

exp(-y). (26) 

where H/^{x) is the forth Hermit polynomial. 



Tip 
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-H^ix) 



B. The slope of the upper critical field in d-wave superconductors 

From Eq.(P^ we obtain that the upper critical field H satisfies 

H° = H - Ar]mde*H^. 

Solving it perturbatively, we get 

H = H° + Ar]mde*{Hy. 

Reinstating the temperature dependence of coefficients of the GL equation, in the simplest 
case a(T) = a'(Tc — T), we finally get 

HiT) = [(Te - T) + 8r]mla' (T, - Tf] . (27) 

We observe that around for a positive rj the H{T) curve bends upwards, in agreement 
with the two field results |T^,[TB|. This effect has been reported in some experiments. However 
one should be cautioned against taking this too seriously. First, as we will discuss in some 
detail later, the coefficient rj should not necessarily be positive in all the samples. For 
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example, twinning is expected to give a negative contribution to it. For negative t] correction 
to the curvature changes sign. Second, ahhough in this study we concentrate on the effects 
of the anisotropy, which is represented by the (scahng) dimension five four derivative term 
d*{Ily — n^)^(i, as we discussed in Section II, there are other rotationally invariant terms of 
the same dimensionality. The second dimension five four derivative term, Td*{lP)'^d, gives 
contribution similar to Eq. (|27|) . One does a calculation similar to that we performed for 
d*{Ily - Ul^d and obtains: 

(T, - T) + 8mla'{r] - r) {T, - Tf . 

This means that for positive r the sign of the correction might be changed. Unlike 77 
which can be fixed by rotation symmetry breaking effects experiments, this rotationally 
invariant correction is more difficult to estimate phenomenologically. 



C. The Abrikosov parameter calculation 

Now we proceed to calculate the Abrikosov's (3a =< > / < MP Here the angular 
bracket is defined as the average over the 2D volume, i.e., < / >= y J (i^r/(r), and V is 
the total volume of the system. If the function /(r) is periodic , it is sufficient to calculate 
the average over one unit cell. 

In the gauge we have chosen in the previous section, a generic solution of the linearized 
equation takes the form \E',fc(x, = exp{iky) ip{x — klfj). Periodicity in y-direction (our 
lattice vector a by assumption is aligned with this axis, see Fig.4) allows the following linear 
combinations: 

n=oo n=oo ^ 2,TTn \ ( I/kI?' \ 

d{x,y)= Cn'^k{x,y)= ^ C„exp(i y\^[x-n ^\ , (28) 

n=-oo n=-oo \ a J \ a ) 

If the second lattice constant is h and it makes an angle a relative to the y axis, the periodicity 
in h direction requires that d{x — 6 sin a,?/ + 6 cos a) = d{x,y) (up to a phase). One can 
achieve it by setting 6 sin a = pAx = p (27rZ|^/a) and Cn+p = Cn exp{i2nnb cos a/a), where p 
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is an integer. For simple Bravais lattices, there is only one vortex in each unit cell. Therefore, 
we can take p = 1. The area of the unit cell is afesina = ^q/H = 27r/|^. 

As a result, all C„ can be fixed up to an overall constant, to be fixed later: 



Cn = exp 
and the wave function becomes: 

ra=oo 

dix, y)= J2 



b nin — 1) 

zTTi — cos a 

a 2 



(29) 



h nin — 1) 

27rz — cos a 

a 2 





27m 


exp 


I y 




a 



il){x — nhsma) 



(30) 



It is convenient to use new rectilinear coordinates whose axes coincide with the vortex 
lattice directions (Fig.4). We shall denote them as X and Y . Their relations to the old x — y 
coordinates are y = Y + X cos a and x = —X sin a. 

The average of \d\'^ is then found by integrating \d\'^ over < X < 6 and < F < a. The 
integration over Y enforces a delta function and simplifies the double summation to 



a6sin 



a sin a / dX\-ilj\( —X ~ nb) sin a]\' 



(31) 



The summation over n converts the integration domain into (— oo, oo). We thus obtain 



b sin a j -oo 
A similar manipulation on {\d\'^) leads to 



dx lihix) 



(32) 



ab sin a 



asina5^+^/,„+„/ exp 



n,rn,n' ,m' 



2Tii - COS a- 
a 



X / dXilj*\{-X - mb) sina] - m'b) sinajipU-X - nb) sina] il)\{-X - n'b) sina] 

Jo 



1 



6 sin a 



m+m' ,n+n' GXp 



n,m,n' ,m' 



2m- COS a 
a 



X 



/ dxtjj*{x — mb sin a) ^|J*{x — m'b sin a) 'i/j{x — nb sin a) il}{x — ra'^sino;) 

J —bsina 



(33) 



Eqs.(|32D and (|33| ) are general expressions for {\d\'^) and We now specialize them 

to our perturbed d field solution Eq.(^). It is easy to see that the correction to starts 
from 77'^. We found that 
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\d\') = 7-^ (l + 
osma V 2 



(34) 



The first order term vanishes, because according to Eg . (|32D , it is proportional to the inner 
product of ipo and ipi. Since we shall be interested only in 0(1]') corrections, we will drop 
this second order term. 

The calculation of {\d\'^) is more involved even in zeroth order Because of the 

presence of the Kronecker delta, there are only three independent summations in Eg. (p^) . 
We choose the summation variables to be 



Z = n + n' = m + m' 
N = n — n' 
M = m — m! 



(35) 



Note that the new discrete variables Z, M and N are not completely independent since they 
have to be either all even or all odd simultaneously. The summation in Eg. (|3^) then becomes 



m+m' ,n+n' 



E E H +11 H H 



(36) 



m.m' ,n,n' even Z even M even A'' odd Z odd A/ odd N 

To zeroth order in rj, the integrand in Eg. (|33|) , after appropriate rearrangement, has a simple 
Gaussian form 



exp 



2sin^a / Z 



''H 



exp 



6^ sin^ a 



My 



N 



(37) 



As before, the summation over Z in Eg.(^) extends the range of the integral over X to 
(— cx), oo), so that the gaussian integral becomes a common factor 



/oo 
dX exp 
-oo 

Pulling out this factor, we obtain 

<\d\'>o=J^^\ E expk.^ 
V 2 sm a , , « 

^^even M,N \ 



2 sin^ a 



/2 



X' 



TT Ih 



2 sin a 



cos a 



exp|27ri- 

odd M,N I ^ 



COS a 



M\2 fN^"^' 



2 y 



6^ sin^ a 



'■H 



(38) 



y) + 



6 sin a 

'■H 



M\2 /A^ 



2 y ^ V 2 
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It is convenient to introduce the complex variable [p5 



C = - exp(ia) = p + i(7. 
a 



(40) 



The coefficients of and A^^ in Eq.(p9[) then become — 27riC*/4 and 2ni(/4:. Finally, we 
obtain 

2' 



Pa = v^{ E exp(27r< 



n 



+ 



E 



(41) 



The above calculation can be straightforwardly extended to include the perturbation of r]. 
The relevant integral now becomes 



dxipi{x — nhsinajipQ {x — msma)ipQ{x — n' sma)ipQ{x — m' sin a) 



where ipi is now given by Eq.(^). The correction of 13a in the first order of 77', after some 
tedious algebra, is: 



7]' \ 

Pa = jv^Re<exp(4i(^) 



exp(27riCn^) (647r^er^rj-^ — 48nan^ + 3) 



/ 1 / / 1 

+ [ n ^ n -\ — ,n n H — 



D. Minimization of the free energy and optimal vortex lattice structure 



Having calculated the Abrikosov parameter I3ai one finds the vortex structure by mini- 
mizing it with respect to (p, p, and a. The minimization with respect to the angle if between 
the vortex lattice and the crystal axes is easily done analytically. The general form of Pa is 



PAi^, p, a) = /?° (p, a)+v \e''^S{p, a) + e-^^^5*(p, a) 



(42) 



Obviously the minimum of (3a is achieved when 



^ = -arg[5(p,cr)]/4±7r/4. 



(43) 



The minimum of (3a is 
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(5T{p,a) = f3\{p,a)-mp,a)\. (44) 

The further minimization of /5™'^(p, cr) is done numerically. In Fig. 5, we show a plot of 
/9™™(p, 0") for 7]' = 0.0193. Due to the fact that the same vortex lattice might be represented 
by several sets of (p, cr), it is enough to consider the region < p < 1/2, see discussion in p5| . 



For every rj', we observed that two degenerate minima appear. One is at p = 1/2, a = 0.663, 
and is clearly an rectangular body centered lattice with a = 53". The corresponding ip 
is zero. Therefore the vortex lattice coincides with the crystal axes. The same result was 
claimed in [0. The other minimum is at p = 0.275 and a = 0.961 and a = 74°, but 
with ip equal to 37°. This corresponds to the previous lattice rotated by 7r/2. To conclude, 
we observed rectangular body centered lattices only. The lowest energy state is doubly 
degenerate. It is interesting to note that the fourfold symmetry of the underlying crystal 
is not completely broken spontaneously by the static vortex lattice - rotations of tt and 
reflections are symmetries of both rectangular lattices. The rj' dependence of a and 
are plotted in Fig. 6 and Fig. 7. We observed that there is a phase transition occured at 
rj' = 0.0235 where the lattice goes continuously from rectangular to square. 

Note that the calculations described in this section, unlike those for the single vortex, are 
valid for arbitrary (not only large k) type II superconductor. Therefore the results can be 
applied to non high Tc materials as well. Using standard methods, one can take into account 
variations of the magnetic field. We do not repeat here the standard relations between the 
free energy and the Abrikosov parameter (3 a P^- One also can calculate corrections to the 



magnetization curves in a standard way using the (3a calculated here. 

Despite the fact that general oblique lattices were considered for the first time for d 
- wave, our numerical analysis shows that they have higher energy than the rectangular 
body centered ones. Intuitively in the symmetric case this is understandable because the 
rectangular lattices are more symmetric. Although for the s - wave superconductors this fact 



has been established , for rotationally non-symmetric superconductors this " argument" is 



not invalid. We are not aware of any mathematical investigations of this question. Moreover, 
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when the vortex lattice starts moving, the rotational symmetry is further explicitly broken. 
As we will see in next section, the general oblique lattices nevertheless are not formed. 



V. THE MOVING LATTICE SOLUTION OF THE TIME DEPENDENT 
GINZBURG-LANDAU EQUATION 

In this chapter we generalize the above procedure to find the structure for a moving 
vortex lattice near ifc2 (the upper critical field itself being a function of the current, as will 
be discussed in Subsection V.B). One can consider the motion as caused by electric field E. 
The vortex lattice velocity is perpendicular to both electric and magnetic field (assumed not 
to be tilted for simplicity and taken to be in the +z direction): E = — v x B. For a general 
direction of the electric field the fourfold symmetry of the system is completely (explicitly) 
broken. Only for several special directions, along the crystal axes [1,0,0], [0,1,0] or along the 
diagonal lines [1,1,0] or [1,1,0] the crystal symmetry is not broken completely, only reduced. 
Even for the simple s - wave time dependent GL equations (TDGL) the problem of finding 
the moving lattice solution is nontrivial. However there exists the "Galilean boost" trick 
^ to solve the linearized (and sometimes even a nonlinear problem for linear response |ll9| ) 



problem. As we will see shortly, for the d - wave equations, even the linearized equation 
does not seem to possess a boosted static solution. 

Technically the steps follow those of the static case. First we find a complete set of 
solutions of the linearized equations using perturbation theory in rj. Then we impose the 
periodicity conditions to construct the vortex lattice wave functions. It is more convenient 
to perform the first step in the gauge aligned in the direction of the electric field, while for 
the second step it is preferable to use a gauge aligned in the direction of the vortex lattice. 
We will combine the two steps using gauge transformation. After the wave function is found, 
it is straightforward to apply the procedure described in the previous section to minimize 
the Abrikosov's [3a- 
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A. The perturbative solution of the hnearized TDGL equations 



To simplify the presentation, we first assume that the direction of the electric field is 
special: along the cryslalline x (or [1,0,0]) direction. In this case the vortices are moving 
in the negative y direction of the coordinate system. We will return to the general case 
afterwards. Now we will construct the perturbative solution to the linearized TDGL Eq.(|TUp 
which we now write in the following "diffusion equation" form: 



(45) 



The vector potential we adopt here is the same as that in section IV, while the electric 
potential can be chosen to be time and y independent: 



-vHx. 



(46) 



In this gauge, the variables t and y trivially separate from x, 

d{x,y,t) = exp{iky) exp{—u!t/'y)4'{x), 

where u can have an imaginary part: u = uor + ioJi. The equation then reduces to a one 
dimensional " Schroedinger type" one (note there is an anti - Hermitian dissipation term) : 

2' 



2mn 



Pi 



^ IT 



Completing the square, rearranging the equation and noting that ui = —ik'yv one obtains: 



2mri 



£ 1 



ri {Ul - Ulf - + ^^'mdvAi;{x) (47) 



(^K - r^y) -ad + -7^mt;^ 



ij{x) = ujRipix). 



(48) 



where a dimensionless quantity g is defined hj g = 'ymdvln and Xq = klfj. The parameter 
ad should be adjusted (i.e., changing the temperature) such that the lowest eigenvalue ur 
becomes zero, otherwise one gets runaway solutions. This is nothing but the Hc2 condition 
generalized to include arbitrary electric field. The operator K defined in Eq.H3) is simply 
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K = Tl^ /2md with x = {x — xo)/Ih shifted by an imaginary amount —ig, so we can write it 



as : 



k = exp{glHPx) Kexp{-g IhPx)- 



(49) 



The perturbation theory to Eq. (|48| ) is most conveniently performed on the shifted field 
defined by: 



il){x) =exp{glHPx)'ip{x) = i){x-ig). 



(50) 



The transformed Hamiltonian is: 



H = exp{-glHPx) (k - r]V^ exp{glHPx) 



Going to the creation and annihilation operators and a representation, Eg. (^81) becomes 

(51) 



a^a H 7]' V{a, a^] 

2 



Here 



V{a, a') = exp 



■ 3 /^t 



a^^ + 1 exp 



■ 3 /^t M 



(52) 



The "potential energy" can be further simplified using the identities: 



exp 
exp 



• 9 M 



a exp 



exp 



■ 9 M 
• 9 /.f M 



a + i 



-\ < ■ 9 
a' + 1—;=. 

V2 



(53) 
(54) 



It is helpful to note that the state resulting from the action of the shifting operator on |0) 
is a coherent state 



-^-^)=exp 



■ 9 /.f M 



|0). 



(55) 



The correction to the eigenvalue ^ (used later to find the phase transition boundary) to the 
first order rj is then easily found: 

1 



^=^-v'i9'-2g' + 2) 



(56) 
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.To the first order in 77, the perturbed ground state is given by: 



4 „/ 



^ = \0) + Y.i-\n){ 



n\ 



n=l 



n 



a' + i—^r + (a + i ^ 



\0)^\0)+v'J2cn\n). (57) 



n=l 



where 



ci = -2V2tg{g^ - 1) 
C2 = -2V2g^ 



C3 = ^ig 
C4 - — 



(58) 



The solution of the original linearized TDGL equation Eq.(^) is simply the shifted iIj{x) 
together with other factors: 



d{x,y,t) = exp[ik{y + vt)] x exp 
1/4 



-^{x kl 



igh 



X 



1 



(59) 



The solution we constructed is restricted to the case when the direction of the electric 
field is along the crystalline x direction. Now we generalize the calculation to arbitrary 
direction of the electric field. In the coordinate system fixed by the vortex lattice (which 
we will use for construction of the vortex lattice solution) the general vortex velocity is: 
Vx = vsm6,Vy = —vcos6, while the electric field is = E cos 6, Ey = EsinO. The angle 
between the crystal [1,0,0] axis and the electric field will be therefore Q = 6 — (f (see Fig. 4). 
The calculation is just a little bit more complicated. It still will be convenient to choose a 
coordinate system in which the direction of the electric field and that of the x axis coincide. 
The perturbed Hamiltonian then becomes: 



V{a, a)) = exp 



• 9 /.f -^ 



(e-2*(^^^)at2 + e2*(^"^)a2)'exp 



■ 9 /^t M 



The corrected solution has the same form as Eq. (|59|) , with the coefficients Cn which now 
depend on the angle {9 — ip), as follows: 
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2 

4^3 



Cl 
C2 
C3 

^ C4 - 2 

The corresponding eigenvalue becomes 



~V2ig [(1 + e-4^(^-^)) _ 2 



^fl + 3e-4^(e-^)) ^2 

4j(0-(^) 



(60) 



ige 



1 



- T] 



COS 4(9 — ip) 4 



2^7^ + 2 



(61) 



B. Dependence of Hc2 on the electric field 

In the simpler case of electric field parallel to one of the crystal axes, from Eg . (1561) the 
new phase boundary equation follows: 

where the second term is a perturbation. All the temperature dependence (at least within the 
confines of the simple Ginzburg - Landau assumption) is contained inside = a'{Tc — T). 
The phase transition line is therefore still quadratic in T, 

Hc2{T) = ho + h{T, - T) + /i2(T, - T)2 (62) 

but first two coefficients have a nontrivial dependance on velocity v: 

2 

h, = 2a'^ (1 - Uvml^h') (63) 

h, = I6a'^r]^ 
e* 

Note that the curvature hasn't changed compared to the static case, but we have two new 
effects. First of all, the electric field (or, equivalently, electric current) shifts Hc2 by a negative 
constant (proportional to i?^ ) to a lower value. This is expected. Secondly, although the 
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curvature /12 doesn't change compared with the static case, the slope hi acquires a negative 
contribution proportional to E"^ . 

In the general case of arbitrary orientation of the electric field, with respect to the crystal 
lattice, only the coefficient ho needs to be changed to: 

i2 



/io = !^ |-1 + [9 + cos 4(e - if)] rim^-f^v^} -fh^ (64) 



We got the interesting result that the shift in Hc2 due to the electric field actually depends 
on the direction of the electic field relative to the crystal lattice. This result should be 
checked experimentally. 

In the s - wave case (or i] = 0) the boundary was first found and discussed in ref. . 
| 20| . There are a couple of peculiarities associated with it like the existence of a metastable 



normal state and the unstable superconductive state. The same applies to the present case. 
As far as we know, these peculiarities haven't been directly observed in low Tc materials. 
It would be interesting to reconsider this question for the high Tc materials. Note also that 
the phase transition is not a usual one (second order which probably turns to weakly first 
order due to fluctuations). In the presence of flux flow the two phases are stationary states 
rather than states in thermal equilibrium. There exist therefore a phase diagram in the 
space containing the current as an external parameter (both magnitude and direction). 



C. Construction of the moving vortex lattice 



Now we would like to follow a procedure similar to that described in Section IV. C for the 
static case to construct a vortex lattice solution. It turns out not to be a straightforward 
generalization. In earlier sections, we used the gauge freedom to make both the scalar and 
the vector potentials independent of y and t. This allows for separation of variables. The 
fact that y variable factored into the form exp{iky) helped us implement the periodicity in 
the y direction (with discrete values of k). However, in general, the vortex lattice will not 
be periodic along this special direction. To construct this general periodic solution, one has 
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to solve a very complicated periodicity constraint equation for the coefficients Ck, where k 
is now a continuous index. 

In the static vortex lattice case, we used the gauge freedom to align the vector potential 
to the vortex lattice. This choice allows us to solve the constraint equation on Ck easily 
since we already had periodicity along the ?/— axis built in. This reduced the problem to a 
discrete one. Furthermore, only a few knS were coupled, and it turned out to be solvable, 
at least for p = 1. This is not the case for the moving vortex lattice. In Subsection A, for 
the problem with electric field and time dependence, we used the gauge freedom to align 
the vector potential with respect to the electric field in order to find the general solution of 
the perturbed Hamiltonian. Now, when we have to use this general solution to construct 
the periodic solution we encounter the problem that we cannot use the gauge freedom to 
simultaneously simplify both problems. Fortunately, in the unperturbed (s-wave) case a 
simple Anzatz for the construction of moving vortex lattice solution exists. This works for 



the linearized TDGL equation with arbitrary direction of the electric field I^Oj. We shall use 
this observation to guide us in obtaining the periodic solution for the moving vortex lattice 
in the presence of perturbation. The solution can be explicitly checked to satisfy TDGL 
equation and the periodicity constraints. 

In subsection A, we were forced to choose an axial gauge as in Eq.(^6l) (which will be 
referred to later as the gauge I), 

A'{x',y',t)=Hx'y' 

y', t)=v A^= -vHx', (65) 

For later convenience, we use x',y' to represent the coordinate in which the electric field is 
along x' direction, while x, y is the coordinate in which the vortex lattice is alligned to the 
y direction (see Fig. 4). The relation is: 

x' = X cos 6 + y sin 6 
' = —X sin 6 + y cos 6 
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Fortunately we can transform our solution to a gauge in which the periodicity is manifest 
and the standard procedure works (referred as gauge II): 



rrid 



A'\x, y,t)^H{x- vj) y + 7— v x z 



= V • = VyH {x - v^t) 



(66) 



where Vx = vsm9,Vy = —vcos9. The gauge transformation between two is determined by 
a phase x(a;, y, i) satisfying 

Vx = - A^^ 
-dtx = - 

One of the solutions is: 



-frudv H 

Y — X H sm ^ cos B 

^ e* 2 



{y' + vtf - x""] + H sin ''9 [x' {y' + vt)] 



(67) 



In this gauge unperturbed lattice can be easily formed using "boosted" solutions, 



1/1 



Vxt knlffj 



(68) 



with standard coefficients: "^^^ = J2n Cn'^i^- These elementary solutions are linearly related 
to the unperturbed normalized eigenfunctions in the gauge I found in Subsection V.A, 



— ' ' exp I j [^^ in' + '^^)] 



(69) 



after the gauge transformation is performed: 

L r°° 



Li /""^ 



Note that gauge transformation and hence quantities Bnk are in general time dependent. 
However, we do not need to keep track of the time dependence here. The reason is that the 
GL equation we are solving with or without the perturbation is time translation invariant. 
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While keeping track of the time dependence of gauge non-invariant quantities Bnk and the 
wave functions is complicated, because in this case we would have to use a time dependent 
gauge choice, the gauge invariant quantities such as Pa are automatically time independent. 
Therefore, to simphfy the calculation, we can set t — 0. 

To find the coefficients Bnf. two gaussian integrations should be performed: 



Bnk = j dx'dy' 



1 



L V ie*^ sin 9 



exp 



1 cos 9 



sin^ 



(70) 
(71) 



2 sine 

We have already found the first order correction to the wave function is gauge I. The 
corresponding expression in the gauge II is: 

L 



" 27T 



where 



S'^l — exp[ik{y' + vt)] x exp 
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dk BnkS'^i 



1 /hV/^ 



H 



X 



X 



( — ) Hn { — - kin - ig 



(72) 



^1 V TT y 

as was shown in the previous section. 

It turns out that after a lengthy calculation, the correction to the wave function in gauge 
II is amazingly simple: 

L roo 



Li /'°° 

ZTT J-oo 



1 / 1 

4 Am6 



1 



H 



X 



knlfl) 



Cm" 



zHr, 



m=l 



X 



- kju - ige 



H 



' A/2^"m! 

An important observation is that the corrected moving lattice solution at t = 

oo 



(73) 



n=— oo 

oo 



1 fH\'/* 



exp {ikny) exp 



■ knl% 



X 



1 + ^7' X! 



m=l 



^2^ 



-.H„ 



X 



- kjH - ige 



-id 



(74) 
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where /c„ = and C„ again given by Eq.(p9D is still invariant under the vortex lattice 
symmetries. From now on we work exclusively in gauge II. Physically this is understood as 
follows. Our perturbation commutes with all the translations generators, in particular with 
the translations defining the lattice structure. Therefore if the unperturbed function has 
certain lattice translation symmetry, the perturbed one will have as well. 



D. The structure and magnetization of the moving lattice 



The standard Abrikosov's procedure to develope an approximation for small order pa- 
rameter around Hc2 can be applied also in the fiux flow case (see pOf). This time however 
the minimization of the Abrikosov parameter (3a does not correspond to minimization of 
energy, but rather to smallest deviation from the exact solution of TDGL equation. The 
"derivation" closely follows the static one. Using the expression for the vortex lattice so- 
lution found in the previous subsection, the correction term in expansion of the Abrikosov 
parameter (3a in 



(3A = (3'A + r]'(3, 



is: 



exp(— 27r2C 



1 , , 1 
2' 2 



Here the function G{n) is defined by 



+ 



(75) 



G{n) = e^'^ieAn^a^n"^ - ASnan^ + 3) 
-Se^^-^c/^ cos2e(87ran2 - 1) 



(76) 



where Q = 6 — ip is the angle between the electric field and the crystal lattice. One immedi- 
ately observes a surprising fact - the dependance on the angle G and velocity v is only via 
the combination g"^ cos 20 where g = '-frridvlH- It factors out as: 
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Pa{v, P, a) = /?o (p, a) + ri Re [e^^'^5(p, a) + cos 2Q^'n\p, a)\ . {11) 

For example, the resulting lattice for O = 7r/4 and arbitrary g will be the same as without 
electric field at all! Also apparent complete breaking of the rotational symmetry by general 
direction of the electric field is not felt by f^A- Indeed, lattice for some arbitrary B and 

is the same as for = and g^' = (7^ cos 20. The fourfold symmetry has been reduced 
however. These results are nontrivial and can be checked experimentally. This degeneracy in 
velocity and in the determination of vortex lattice may be a refiection of some dynamical 
symmetries which we have so far failed to see yet. 

We get the e^^**^ harmonics in Eq. ( [77|) in addition to the fourth harmonic that appeared 

in the static case. The minimization with respect to (f still can be done analytically, although 

the algebraic equation in this case is quartic. For fixed rj' and (7^ cos 20, the minimization 

with respect to p and a was performed numerically and we again obtain only rectangular 

body centered lattices aligned either to the crystalline axes. The angle a turns out to be 

only weakly dependent on the combination g"^ cos 20. For example, for positive f]' = 0.015 , 

we obtain a = a{g = 0) + 1.0(7^1 cos20| (in degrees) where a{g = 0) = 69.3°. The Abrikosov 

13a is simply related to the slope of the magnetization curve 

. dM 1 
-"'dif = (2.^ - 1)0, 

as well as to other thermodynamic quantities. All of them therefore exhibit very simple 
dependence on the velocity v , or, equivalently, on the current J. 

The fact that the optimal lattice is rectangular body centered is a bit mysterious. Rota- 
tional symmetry is completely broken by both the electric field and by the underlying crystal 
lattice. It is not easy to attribute the advantage of this lattice structure to some simple phys- 
ical origin. It might be that it is a consequence of using the Abrikosov approximation, and 
therefore beyond this approximation the lattices might not be rectangular. Note also that 
it was also surprising that (3a was independent of the orientation of the electric field even 
in the s - wave calculation |2^. As far as we know, the preferred orientation of the moving 



lattice has not been observed in either low Tc or high Tc type II superconductors. 
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VI. THE NON-LINEAR CONDUCTIVITY NEAR Hc2 



In this section we consider the dissipation in vortex cores due to flux flow. As it is well 
known, the fourfold symmetry forces the conductivity tensor cijj, defined by J, = CijEj, to be 
rotationally symmetric. Namely, aij = a5ij+a^ Sij. Here a is the usual (Ohmic) conductivity, 

is the Hall conductivity and Sij is the antisymmetric tensor. The additional term in the 
free energy corrects the values of a and cr^, but the correction is of the order rj and therefore 
small. So, to see anisotropy, we definitely would like to go beyond linear response. This has 
been done for simple s - wave TDGL ||20| near Hc2 ■ We will neglect pinning and consider 
motion of a very large bundle. While there is a normal component of the conductivity, here 
we will concentrate on the contribution of the supercurrent only. For the discussion of the 



relative contribution of the two see 20 



A. Condensate for the moving lattice 

To calculate the transport properties due to flux flow, we have to compute two quan- 
tities. The first is the expression for the electric current and will be obtained in the next 
subsection. The second is the normalization factor for the d - wave order parameter. In the 
previous sections we only needed Abrikosov's 13a which is insensitive to the overall scale of 
the condensate, but now we need to calculate it. 

Here \I'„ and 5\E'„ have been calculated in the previous section, is the normalization 
and A^o is the normalization of . The calculation is standard. We again expand it to first 
order in t]'. The normalization is determined from the minimization of the free energy as 

< d*d >= 

2(3 (3a 

where < ... > denote the space average, ad and (3 are coefficients of the GL equation. The 
Abrikosov parameter (3a has its own rj' expansion: (3a = (3\ + v'Pa calculated in Section 
VI. C. Combining the two one obtains 
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1 



which will be used to calculate the current. 



1 / 



+ 



2Re < >' 



B. The direct and the Hall currents 



The anisotropy term in the free energy Eq.(^ contains four covariant derivatives and 
consequently the electric current, in addition to the usual expression, contains additional 
terms. The leading order current is given by 

■{d*{nd) + {nd)*d) 



2m, 



The anisotropic perturbation to the current is 



J\d) = e*ri±" 



-e r]y 



n 



//2 TT"2 



u':']d 



iKd) + [K (nf -nf)rf]*rf + c.c. 



{K - nf ) d]* (n;'d) + [n'; (n'^^ - nf ) d 



d + c.c. 



(79) 



where 11'^ = cosipUx + sim^IIj^^ Uy = cos ipUy — siny^IIj;^ and x" = cos (/^x + sin y^y, y" = 
cosy^y — sini^x (See Fig. 4). 

Substituting the condensate d = N {"^ + ?7'5\E') with \E' and 5\E' determined in the previous 
subsection, one obtains the expansion of J" to the first order in rj': 

j^^N^ (2^) {[(^*n^) + (^n^*)] + 

where Sji comes from the correction to A^^, and 6^2 contains the correction to the wave 
function 6"^. 



(80) 



5h 
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_ / e*ad 1 \ 4Re(5**n*) 
The expansion to first order in 77' of J is: 



(81) 



e*ad 1 



'3l 



X 



2/5 

;(n;;^ - nf ) (n^^) + [n^ (n;;^ - nf ) ^ 



^ + c.c. 



-y 



< > 



Summing up the corrections, one obtains the corrected current as : 

J = J« + J'' = j + r7'<5j = j + V (5ji + 5h + 5h) 

Here we have used the fact that the operator 11 is Hermitian, and the curly bracket denotes 
anticommutator. Note that although the total wave function is a linear combination of 
^n{x-i y), after averaging over the 2D space all the components decouple due to the exp{ikny) 
factors and the fact that the current is quadratic in This means that it suffices to consider 
only one of the components to calculate J. The results for the three contributions to the 
current correction are: 



-2 Re < > e*ad 1 

~ J — 



< > 



P PI 



-9 



e*ad 1 4 Re < 5**n* > e*ad 1 



4/3md < > 



P P% 



1 + - cos 46 



1 + - cos 46 



+ 4/ 1 (7v X z) 
- - 2| (7v X z) 



(82) 



(83) 



e*adil 1 1^,2 Re (v[/*{n;'2-nf,n^_^^^ 21^x1/* {n;;2-nf,n;'}v[/ 



2/3 rudPA \ < > " < > 

= { (1 + + 2I (7V X z) - / sin 46 (7V) } 

P Pa 

where Q = 9 — if as before. Summing them up we got our final expression: 



(84) 



33 



5j = -§j + '^^9' (1 + cos4e) (7v X z) 

HA P Pa 

e*ad 1 2 • 



sin4e(7v) 



/3 P'a 

From this, one obtains the simple resuhs Eq.(^^ advertised earher. Note that all the 
terms are cancelled. 



VII. CONCLUDING REMARKS 

Instead of summarizing the results (which has been done in Section I), we briefly comment 
on the possibility of observation of various phenomena quantitatively discussed in this paper. 

.1. Internal structure of a single anisotropic vortex Although direct observation 

of the order parameter using scanning - tunneling- microscopy (STM) or the magnetic 
field distribution using electron holography |^ or other techniques is possible, the de- 
tailed effects hotly debated by theoreticians (where the zeroes of the s field are located, 
small distance asymptotics) probably do not have a significant impact on such experiments. 
One also should note that the Ginzburg - Landau framework adopted here might not be 
applicable close to the vortex center where microscopic excitation spectrum becomes impor- 
tant. An approach using elements of the microscopic theory (via Bogoliubov - DeGennes 



equations along the lines of and ||T5[) will be necessary. 



.2. Structure of static and moving anisotropic vortex lattice Static vortex 
lattice has been observed using small angle neutron scattering ||^ and tunneling spectroscopy 
0. Although moving vortices have been directly observed using electron tomography |27], to 



our knowledge the shape and orientation of moving large bundles hasn't been observed as 
yet. Moving vortex lattice is much more sensitive to pinning effects than the static lattice. 
In the static case pinning can just slightly distort or cause breakup of the crystal to smaller 
peaces. For the moving flux lattice pinning is expected to be much more significant. The 
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orientation effect that we predict is very small, but the asymmetry in magnetization might 
be quite significant. 

We found the transtion point for the parameter rj = 7]mde*H is at r/^ = 0.0235. This 
transition between the rectangular and the square lattices might be seen in neutron scattering 
experiments, since the square lattice has higher symmetry (number of spots is reduced to 
four at the transition). Note that by increasing magnetic field the critical rj'^ can be exceeded 
without changing the sample (77 is independent of magnetic field). 

.3. Static transition to the normal phase As is well known, in the presence 
of fluctuations, the second order phase transition from superconducting to normal state 
becomes a weakly first order melting line into the vortex liquid. This is the reason that the 
present study of the diagram will be useless for BSCCO which has a relatively large Ginzburg 
number. For YBCO and the low temperature superconductors, the curvature of the phase 
transition line can in principle provide an estimate of 1] with reservations mentioned in the 
end of Subsection IV. C. 

.4. Dynamic phase diagram Dynamical phase diagram, namely transition from the 
moving lattice to normal (or moving liquid) should be complicated by pinning effects. How- 
ever, provided these could be overcome a number of interesting effects could be observed. 
First, even neglecting the rotational symmetry breaking effects, there are a number of pecu- 
liarities associated with normal - superconducting boundary noted by Thompson and Hu like 
the existence of a metastable normal state and the unstable superconductive state. As far 
as we know, these peculiarities haven't been convincingly observed in low materials. 
It would be interesting to reconsider this question for the high Tc materials. 

In the s - wave case however, the phase diagram cannot depend on orientation of the 
current. We calculated this orientation dependance on the angle between the atomic lattice 
and the direction of current or electric field to first order in t], see Eq. (|6^j63|j6^) . New effects 
include the change in slope of Hc2 as function of temperature, not only in curvature. 
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a. Nonlinear I-V curves and magnetization One should be able to measure currents in 
the same sample oriented differently with respect to the atomic crystal. Note that the effect 
can be seen in low temperature anisotropic superconductors, not necessarily in YBCO. The 
simplicity of the expressions for both direct and Hall current Eqs. calls for some special 
ways to verify it experimentally. The angular dependence of the magnetization near the 
transition (given by Eqs. (|78| , [75|) ) might be large enough to be measurable. 

There are number of limitations of our approach which can be lifted by possible exten- 
sions. One of them is the assumption of exact fourfold symmetry. Deviations from it in a 
form of different coefficients of the gradient terms in x and y directions have already been 
studied recently |^ using the two field formalism. If they happen to be small, they can 
be easily added perturbatively. These effects of explicit breaking are clearly quite different 
from those of spontaneous breaking of the fourfold symmetry studied here. Our results for 
the latices are limited to fields close to Hc2 only. It is possible, although more difficult to 
extend them to lower magnetic fields. Another interesting direction is the influence of the 
anisotropy on vortex fluctuations in the lattice. We hope to address these issues in the 
future. 

Also anisotropy influences fluctuations not considered here. In addition, the effective one 
component approach allows to consider possibilities not apparent within the two field one. 
For example, the coefficient t], in principle, can be negative despite the fact that within the 
two field formalism it should be positive. Twinning is expected to reduce the value of the 
parameter. 

APPENDIX: COMPARISON OF THE TWO FIELD AND THE ONE FIELD 

RESULTS FOR SINGLE VORTEX 

In the two field formulation, the small r asymptotics of the solution for the d - wave 
component of the isolated vortex is given by : 

d{r, 0) ^ {dir + c/gr^) e''^ (Al) 
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where the subleading term coefficients is 



= - 



A. 



1 + 



ho 



di 



neglecting terms proportional to ho/Hc2{0) ||T3|. The s-component asymptotics is 



sir, 



--(n;-n^Kr) 



dire 



-i<f> 



(A2) 



2 V«^d/ 

The expression Eq.(^^) is different from what was obtained in Eq. (22), but qualita- 
tively the behavior is not affected. We also found that it doesn't follow from their Eq. (19), 
because to the same order of approximation had a non-vanishing term proportional to 
g3i0_ Nevertheless the concluding statement in [|13[ is basically correct. Following the same 
argument leading to an estimate of the maximal amplitude of s(r) as in we obtained 



Iv 



do 4 \asQt, 



apparently this correction accounts for the 20% error cited in |T^ . 

The asymptotic form of the wave function was used to make a topological argument 
about poles in the s - wave. Due to the different winding number of small r and large r 
asymptotics of s(r, 0), there must exist four poles in the intermediate region. This was shown 
numerically in [T^. Ting et al [T^ however performed a similar calculation, but didn't get 
the poles. Our calculation, which is much simpler than the two field one, conffims the former 
and shows clearly four poles on the x and y axis, independent of what kind of approximate 



(i— component wave function one chooses. We suspect that the numerical simulation in [17 



was not sensitive enough to resolve these poles |]36[| 
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Figure Captions 



FIG. 1 A single vortex solution of the one component GL equation. The coefficient function 
/i) /-3) /5 for the harmonics e*"^, e"^*"^, e^*"^, respectively, /j's are given in units of \E'o = ad/2/3, 
and r is given in units of ^d- See Eq. (|TB|),(|T^),(|I^) 

FIG. 2 The d-field of a single vortex for rj = 0.15. Only the absolute value of the d field in 
units of \Efo is shown, (a) Contour plot, (b) Three dimensional plot. 

FIG. 3 The s-field of a single vortex for t] = 0.15. (a) Contour plot, (b) Three dimensional 
plot. Note that there are four singularities on which the s-field vanishes. 

Fig. 4 The coordinate system used in our calculations, this defines the angles 6, (f and 0. 

FIG. 5 The Abrikosov parameter (3^ as a function of the lattice parameters (p, a) (defined 
in Eq.(^Dp). There are three degenerate local minima. Oblique lattices are on the lines 
p = 1/2 and p^ + o"^ = 1. The two points A and B are related by p 1/p and therefore 
represents the same lattice. Point C represents the same rectangular lattice rotated by 90°. 

FIG. 6 The angle a as a function of rj', the two branches correspond to lattices related by 
a rotation of 90°. A continuous transition from the rectangular lattice to the square lattice 
happens at r]'^ = 0.235. 

FIG. 7 The Abrikosov parameter jS^ as a function of r]' for triangular, square and optimal 
rectangular body centered lattices,respectively. At the transition point rj'^, the rectangular 
lattice is taken over by the square lattice. Note that rj' is proportional to the magenetic field 
H. 
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